# ===============================================================#
#                     Replication files for:                     #
#.  "Attitudinal and Behavioral Legacies of Wartime Violence:    #
#                      A Meta-Analysis"                          #
#                        Joan Barceló                            #
#               American Political Science Review                #
#               Last update: September 3, 2025                   #
# ===============================================================#

#################################
# Figure 3: Figure on Outgroup Hostility and Ingroup Favoritism
################################

## ------------------------ Load Packages ------------------------

library(tidyverse)
library(ggplot2)
library(patchwork)
library(grid)

# ---------- Paths ----------

base_path_meta    <- "~/Datasets/"
base_path_effects <- "~/Results/"

read_if <- function(path) if (file.exists(path)) read.csv(path) else NULL
lbl <- function(b, p) paste0(round(b, 2), ifelse(p < .001, "***", ifelse(p < .01, "**", ifelse(p < .05, "*", ""))))

# Styles
col_full <- "indianred3"
col_biv  <- "steelblue4"
col_qe   <- "forestgreen"
sh_full  <- 15
sh_biv   <- 15
sh_qe    <- 17

## ------------------------ Load datasets ------------------------
# Outgroup hostility, three outcomes
meta.warenemies        <- read_if(paste0(base_path_meta, "meta.warenemies.csv"))
meta.warenemies.random <- read_if(paste0(base_path_meta, "meta.warenemies_random.csv"))
meta.warenemies_biv    <- read_if(paste0(base_path_meta, "meta.warenemies_biv.csv"))
warenemies.effects     <- read_if(paste0(base_path_effects, "9.results_wartime_enemies_models.csv"))
if (is.null(warenemies.effects)) warenemies.effects <- read_if(paste0(base_path_effects, "9.results_warenemies_models.csv"))

meta.threat            <- read_if(paste0(base_path_meta, "meta.threat.csv"))
meta.threat.random     <- read_if(paste0(base_path_meta, "meta.threat_random.csv"))
meta.threat_biv        <- read_if(paste0(base_path_meta, "meta.threat_biv.csv"))
threat.effects         <- read_if(paste0(base_path_effects, "11.results_threat_models.csv"))

meta.intergroup        <- read_if(paste0(base_path_meta, "meta.intergroup.csv"))
meta.intergroup.random <- read_if(paste0(base_path_meta, "meta.intergroup_random.csv"))
meta.intergroup_biv    <- read_if(paste0(base_path_meta, "meta.intergroup_biv.csv"))
intergroup.effects     <- read_if(paste0(base_path_effects, "10.results_intergroup_models.csv"))

# Ingroup identification and favoritism, three outcomes
meta.groupvoting        <- read_if(paste0(base_path_meta, "meta.groupvoting.csv"))
meta.groupvoting.random <- read_if(paste0(base_path_meta, "meta.groupvoting_random.csv"))
meta.groupvoting_biv    <- read_if(paste0(base_path_meta, "meta.groupvoting_biv.csv"))
groupvoting.effects     <- read_if(paste0(base_path_effects, "13.results_groupvoting_models.csv"))

meta.groupid        <- read_if(paste0(base_path_meta, "meta.groupid.csv"))
meta.groupid.random <- read_if(paste0(base_path_meta, "meta.groupid_random.csv"))
meta.groupid_biv    <- read_if(paste0(base_path_meta, "meta.groupid_biv.csv"))
groupid.effects     <- read_if(paste0(base_path_effects, "14.results_groupid_models.csv"))

meta.ingroup_trust        <- read_if(paste0(base_path_meta, "meta.ingroup_trust.csv"))
meta.ingroup_trust.random <- read_if(paste0(base_path_meta, "meta.ingroup_trust_random.csv"))
meta.ingroup_trust_biv    <- read_if(paste0(base_path_meta, "meta.ingroup_trust_biv.csv"))
ingroup_trust.effects     <- read_if(paste0(base_path_effects, "12.results_ingroup_trust_models.csv"))

# Stack effects
re3_estimates <- bind_rows(
  warenemies.effects, threat.effects, intergroup.effects,
  groupvoting.effects, groupid.effects, ingroup_trust.effects
)

## ======================== PANEL A: Attitudinal hardening ========================

## A1) Antagonism toward wartime enemies
origN <- n_distinct(meta.warenemies$authoryear); origK <- nrow(meta.warenemies)
qeN   <- if (is.null(meta.warenemies.random)) 0 else n_distinct(meta.warenemies.random$authoryear)
qeK   <- if (is.null(meta.warenemies.random)) 0 else nrow(meta.warenemies.random)

pd_warenemies <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Antagonism toward wartime enemies")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model %in% c("Full - All", "Original - All") ~ 1.012,
      sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental") ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model %in% c("Full - All", "Original - All") ~ col_full,
      sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental") ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental"), sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% tidyr::drop_na(y)

warenemies.plot <- ggplot(meta.warenemies, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.warenemies_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .2, width = 0, height = .012) +
  geom_point(data = pd_warenemies, aes(x = coef, y = y, color = col), shape = pd_warenemies$shp, size = 2.6, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = pd_warenemies, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .9) +
  geom_text(data = pd_warenemies, aes(x = coef, y = y + .005, label = lab, color = col), size = 3, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "", y = "Antagonism\ntoward\nwartime enemies") +
  theme_void() +
  theme(axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
        plot.margin = margin(0,0,0,0), legend.position = "none") +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", origN, "\nk = ", origK), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qeN,   "\nk = ", qeK ), size = 3)

## A2) Perceived threat from wartime enemy
origN <- n_distinct(meta.threat$authoryear); origK <- nrow(meta.threat)
qeN   <- if (is.null(meta.threat.random)) 0 else n_distinct(meta.threat.random$authoryear)
qeK   <- if (is.null(meta.threat.random)) 0 else nrow(meta.threat.random)

pd_threat <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Perceived threat from wartime enemy")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model %in% c("Full - All", "Original - All") ~ 1.012,
      sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental") ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model %in% c("Full - All", "Original - All") ~ col_full,
      sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental") ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental"), sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% tidyr::drop_na(y)

threat.plot <- ggplot(meta.threat, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.threat_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .2, width = 0, height = .012) +
  geom_point(data = pd_threat, aes(x = coef, y = y, color = col), shape = pd_threat$shp, size = 2.6, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = pd_threat, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .9) +
  geom_text(data = pd_threat, aes(x = coef, y = y + .005, label = lab, color = col), size = 3, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "", y = "Threat\nperceptions") +
  theme_void() +
  theme(axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
        plot.margin = margin(0,0,0,0), legend.position = "none") +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", origN, "\nk = ", origK), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qeN,   "\nk = ", qeK ), size = 3)

## A3) Intergroup distrust, bias or discrimination
origN <- n_distinct(meta.intergroup$authoryear); origK <- nrow(meta.intergroup)
qeN   <- if (is.null(meta.intergroup.random)) 0 else n_distinct(meta.intergroup.random$authoryear)
qeK   <- if (is.null(meta.intergroup.random)) 0 else nrow(meta.intergroup.random)

pd_intergroup <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Intergroup distrust, bias, or discrimination")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model %in% c("Full - All", "Original - All") ~ 1.012,
      sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental") ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model %in% c("Full - All", "Original - All") ~ col_full,
      sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental") ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental"), sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% tidyr::drop_na(y)

intergroup.plot <- ggplot(meta.intergroup, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.intergroup_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .2, width = 0, height = .012) +
  geom_point(data = pd_intergroup, aes(x = coef, y = y, color = col), shape = pd_intergroup$shp, size = 2.6, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = pd_intergroup, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .9) +
  geom_text(data = pd_intergroup, aes(x = coef, y = y + .005, label = lab, color = col), size = 3, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "Effect size", y = "Intergroup\ndistrust,\nbias, or\ndiscrimination") +
  theme_minimal() +
  theme(
    axis.title.x = element_text(size = 10, face = "bold"),
    axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
    axis.text.y = element_blank(), axis.ticks.y = element_blank(),
    panel.grid = element_blank(),
    plot.margin = margin(0,0,0,0), legend.position = "none"
  ) +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", origN, "\nk = ", origK), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qeN,   "\nk = ", qeK ), size = 3)

## ======================== PANEL B: Ingroup identification and favoritism ========================

## B1) Group-based Voting
origN <- n_distinct(meta.groupvoting$authoryear); origK <- nrow(meta.groupvoting)
qeN   <- if (is.null(meta.groupvoting.random)) 0 else n_distinct(meta.groupvoting.random$authoryear)
qeK   <- if (is.null(meta.groupvoting.random)) 0 else nrow(meta.groupvoting.random)

pd_groupvoting <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Group Voting")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model %in% c("Full - All", "Original - All") ~ 1.012,
      sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental") ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model %in% c("Full - All", "Original - All") ~ col_full,
      sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental") ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental"), sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% tidyr::drop_na(y)

groupvoting.plot <- ggplot(meta.groupvoting, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.groupvoting_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .2, width = 0, height = .012) +
  geom_point(data = pd_groupvoting, aes(x = coef, y = y, color = col), shape = pd_groupvoting$shp, size = 2.6, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = pd_groupvoting, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .9) +
  geom_text(data = pd_groupvoting, aes(x = coef, y = y + .005, label = lab, color = col), size = 3, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "", y = "Group-based\nvoting") +
  theme_void() +
  theme(axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
        plot.margin = margin(0,0,0,0), legend.position = "none") +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", origN, "\nk = ", origK), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qeN,   "\nk = ", qeK ), size = 3)

## B2) Group Identification
origN <- n_distinct(meta.groupid$authoryear); origK <- nrow(meta.groupid)
qeN   <- if (is.null(meta.groupid.random)) 0 else n_distinct(meta.groupid.random$authoryear)
qeK   <- if (is.null(meta.groupid.random)) 0 else nrow(meta.groupid.random)

pd_groupid <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Group Identification")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model %in% c("Full - All", "Original - All") ~ 1.012,
      sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental") ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model %in% c("Full - All", "Original - All") ~ col_full,
      sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental") ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental"), sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% tidyr::drop_na(y)

groupid.plot <- ggplot(meta.groupid, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.groupid_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .2, width = 0, height = .012) +
  geom_point(data = pd_groupid, aes(x = coef, y = y, color = col), shape = pd_groupid$shp, size = 2.6, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = pd_groupid, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .9) +
  geom_text(data = pd_groupid, aes(x = coef, y = y + .005, label = lab, color = col), size = 3, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "", y = "Group\nidentification") +
  theme_void() +
  theme(axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
        plot.margin = margin(0,0,0,0), legend.position = "none") +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", origN, "\nk = ", origK), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qeN,   "\nk = ", qeK ), size = 3)

## B3) Ingroup trust
origN <- n_distinct(meta.ingroup_trust$authoryear); origK <- nrow(meta.ingroup_trust)
qeN   <- if (is.null(meta.ingroup_trust.random)) 0 else n_distinct(meta.ingroup_trust.random$authoryear)
qeK   <- if (is.null(meta.ingroup_trust.random)) 0 else nrow(meta.ingroup_trust.random)

pd_ingroup <- re3_estimates %>%
  filter(tolower(outcome) == tolower("Ingroup trust")) %>%
  transmute(
    sample_model = paste(model_type, sample_type, sep = " - "),
    coef, ci.lb95, ci.ub95, pval,
    y = case_when(
      sample_model == "Bivariate - All" ~ 1.000,
      sample_model %in% c("Full - All", "Original - All") ~ 1.012,
      sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental") ~ 0.988,
      TRUE ~ NA_real_
    ),
    col = case_when(
      sample_model == "Bivariate - All" ~ col_biv,
      sample_model %in% c("Full - All", "Original - All") ~ col_full,
      sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental") ~ col_qe,
      TRUE ~ NA_character_
    ),
    shp = ifelse(sample_model %in% c("Full - Quasi-experimental", "Original - Quasi-experimental"), sh_qe, sh_full),
    lab = lbl(coef, pval)
  ) %>% tidyr::drop_na(y)

ingroup_trust.plot <- ggplot(meta.ingroup_trust, aes(x = coef, y = 1, size = 1/se)) +
  geom_jitter(shape = 21, fill = "gray70", color = "gray70", alpha = .3, width = 0, height = .012) +
  geom_jitter(data = meta.ingroup_trust_biv, aes(x = coef, y = 1, size = 1/se),
              shape = 4, stroke = 1, color = "gray70", alpha = .2, width = 0, height = .012) +
  geom_point(data = pd_ingroup, aes(x = coef, y = y, color = col), shape = pd_ingroup$shp, size = 2.6, stroke = 1.5, fill = "white") +
  geom_errorbarh(data = pd_ingroup, aes(y = y, xmin = ci.lb95, xmax = ci.ub95, color = col), height = 0, size = .9) +
  geom_text(data = pd_ingroup, aes(x = coef, y = y + .005, label = lab, color = col), size = 3, fontface = "bold") +
  scale_color_identity() + scale_size_continuous(range = c(1, 10), guide = "none") +
  labs(x = "Effect size", y = "Ingroup\ntrust") +
  theme_minimal() +
  theme(
    axis.title.x = element_text(size = 10, face = "bold"),
    axis.title.y = element_text(angle = 90, size = 10, face = "bold"),
    axis.text.y = element_blank(), axis.ticks.y = element_blank(),
    panel.grid = element_blank(),
    plot.margin = margin(0,0,0,0), legend.position = "none"
  ) +
  xlim(c(-.4, .4)) + ylim(c(.985, 1.025)) +
  ggplot2::annotate("text", x = .30, y = 1.000, label = paste0("All:\nN = ", origN, "\nk = ", origK), size = 3) +
  ggplot2::annotate("text", x = .40, y = 1.000, label = paste0("QE:\nN = ", qeN,   "\nk = ", qeK ), size = 3)

## ------------------------ Legend and titles ------------------------
legend_plot <- ggplot() +
  xlim(0, 1) + ylim(0, 1) + theme_void() +
  ggplot2::annotate("point", x = 0.70, y = 0.90, shape = sh_full, size = 3, color = col_full, fill = col_full) +
  ggplot2::annotate("text",  x = 0.73, y = 0.90, label = "Full models", hjust = 0, vjust = 0.5, size = 3) +
  ggplot2::annotate("point", x = 0.70, y = 0.70, shape = sh_biv,  size = 3, color = col_biv,  fill = col_biv) +
  ggplot2::annotate("text",  x = 0.73, y = 0.70, label = "Bivariate models", hjust = 0, vjust = 0.5, size = 3) +
  ggplot2::annotate("point", x = 0.70, y = 0.50, shape = sh_qe,   size = 3, color = col_qe,  fill = col_qe) +
  ggplot2::annotate("text",  x = 0.73, y = 0.50, label = "Quasi experimental models", hjust = 0, vjust = 0.5, size = 3) +
  theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  labs(x = "", y = "")

panel_title <- function(title) {
  patchwork::wrap_elements(grid::textGrob(title, gp = gpar(fontsize = 12, fontface = "bold"), just = "center"))
}

## ------------------------ Compose and export Figure ------------------------
Figure3 <- (
  panel_title("Panel A. Attitudinal Hardening Toward Wartime Enemies") /
    warenemies.plot /
    threat.plot /
    intergroup.plot /
    panel_title("Panel B. Ingroup Identification and Favoritism") /
    groupvoting.plot /
    groupid.plot /
    ingroup_trust.plot /
    legend_plot
) & theme(axis.title.x = element_text(face = "bold")) & plot_layout(guides = "collect")

#out_path <- "/Users/jbs548/Library/CloudStorage/Dropbox/Joan Barcelo/Present/NYUAD Assistant Professor/Research/Papers/Work in Progress/Meta_analysis/Analysis/Replication_materials_APSR/Figures/"
#pdf(paste0(out_path, "Figure3.pdf"), height = 15, width = 8.4, paper = "letter")
#print(Figure3)
#grid.draw(linesGrob(x = unit(c(0.525, 0.525), "npc"),
#                    y = unit(c(0.60, 0.92), "npc"),
#                    gp = gpar(col = "darkgray", alpha = 0.4)))
#grid.draw(linesGrob(x = unit(c(0.525, 0.525), "npc"),
#                    y = unit(c(0.18, 0.50), "npc"),
#                    gp = gpar(col = "darkgray", alpha = 0.4)))
#dev.off()

